Przed oddaniem zadania upewnij się, że wszystko działa poprawnie. Uruchom ponownie kernel (z paska menu: Kernel$\rightarrow$Restart) a następnie wykonaj wszystkie komórki (z paska menu: Cell$\rightarrow$Run All).
Upewnij się, że wypełniłeś wszystkie pola TU WPISZ KOD lub TU WPISZ ODPOWIEDŹ, oraz
że podałeś swoje imię i nazwisko poniżej:
NAME = "MICHAL MARCINIAK"
import pandas as pd
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn import linear_model
import pickle
import numpy as np
from matplotlib import pyplot as plt
plt.rcParams["animation.html"] = "jshtml"
%matplotlib inline
from pprint import pprint
with open('./data/fish-preprocessed.pkl', 'rb') as f:
data = pickle.load(f)
train_data = data['train']
test_data = data['test']
Podejść do "wektoryzacji" zmiennych kategorycznych jest wiele i ich stosowanie jest niezbędne w przypadku większości modeli. W przypadku modeli liniowych dobór odpowiedniego podejścia jest jednak wyjątkowo istotny, ponieważ model nie odtworzy ukrytych zależności - w przeciwieństwie do np. sieci neuronowych (przynajmniej w teorii). Oto wybór możliwości:
Osobny model dla każdej z kategorii
Mało wydajne, ponieważ zwiększamy liczbę stopni swobody (VC-dimensiom, klątwa wymiarowości itp) przy tej samej liczbie danych. Ale jeżeli każda z klas posiada inny proces generujący, to po co łączyć w jeden model?
Utworzenie $k$ bądź $k - 1$ nowych zmiennych niezależnych (gdzie $k$ to liczba kategorii)
Popularne i łatwe w użyciu podejście, nie wymaga zrozumienia zbioru danych, daje dobre wyniki dla modeli nie-liniowych. Można przeprowadzić na wiele sposobów - https://stats.idre.ucla.edu/spss/faq/coding-systems-for-categorical-variables-in-regression-analysis-2/
Warunkowanie istniejących zmiennych niezależnych
Wymaga zrozumienia zbioru danych (dla jakich zmiennych niezależnych wpływ na zmienną zależną może być różny zależnie od kategorii), ale dla modeli liniowych będzie dawać najlepsze wyniki (jeżeli przeprowadzone prawidłowo = zgodnie z procesem generatywnym/z jego najlepszą liniową aproksymacją).
Poza trudnością w identyfikacji warunkowania, pojawiają się również trudności implementacyjne. Ręczne przygotowanie zbioru danych w ten sposób jest czasochłonne i łatwo o błąd. Bardzo dobrym rozwiązaniem jest podejście z pakietu R, tzw. formula function, którą posiada bardzo funkcjonalną składnię. W Pythonie z tej samej składni można skorzystać w pakiecie statsmodels.
Zastosujemy rozwiązanie nr 2. w wersji dummy / one-hot / one-of-K coding, czyli każda z klas będziała miała swoją kolumnę wypełnioną zerami i jedynkami. Dla jednej z klas (klasy referencyjnej) moglibyśmy nie tworzyć kolumny, ponieważ wyraz wolny "przejąłby" referencyjny wpływ na zmienną zależną. Gdybyśmy stosowali model inny niż liniowy z wyrazem wolnym takie podejście mogłoby się nie sprawdzić. Wyraz wolny, wkrtóce zniknie z naszego modelu liniowego, więc pozostawiamy $k$ kolumn.
one_hot = OneHotEncoder(drop=None, sparse=False)
train_data = train_data.join(pd.DataFrame(
one_hot.fit_transform(train_data['Species'].to_numpy().reshape(-1, 1)),
columns=[f'species_{i}' for i in range(7)],
index=train_data.index,
))
train_data
A teraz, wcześnie dopasowany na danych treningowych transformator należy zaaplikować na danych testowych:
test_data = test_data.join(pd.DataFrame(
one_hot.transform(test_data['Species'].to_numpy().reshape(-1, 1)),
columns=[f'species_{i}' for i in range(7)],
index=test_data.index,
))
Oryginalną kolumnę zostawiamy w DataFrameach, ponieważ ułatwi nam to wizualizację.
Model: $$ y_i = \pmb{x}_i \pmb{\beta} $$
Regularyzacja $L_2$:
$$ \mathcal{R}(\pmb{\beta}) = \frac{\lambda}{2} \left\|\pmb{\beta}\right\|_2^2,$$gdzie $\lambda$ to parametr metody regulujący siłę kary za złożoność.
Funkcja kosztu: $$ S(\pmb{\beta}) = \frac{1}{2} \sum_{i=1}^N \left| y_i - (\beta_0 + \pmb{x}_i \pmb{\beta}) \right|^2 + \frac{\lambda}{2} \left\|\pmb{\beta}\right\|^2, $$ gdzie macierz $X$ nie posiada tutaj dodatkowej kolumny $1$ odpowiadającej za wyraz wolny. Jest on traktowany osobno, ponieważ nie chcemy na niego zakładać regularyzacji.
Estymacja:
Analityczna Najpierw musimy dopasować wyraz wolny $\beta_0$. Ma być to wartość o najmniejszym błędzie średniokwadratowym, gdyby wszystkie cechy miały wartość 0.
$$ \hat{\beta}_0 = \frac{1}{N} \sum_{i=1}^N y_i $$
Następnie wyraz wolny należy odjąć od zmiennej niezależnej co daje następującą funkcję kosztu: $$ S(\pmb{\beta}) = \frac{1}{2} \sum_{i=1}^N \left| (y_i - \hat{\beta}_0) - \pmb{x}_i \pmb{\beta}) \right|^2 + \frac{\lambda}{2} \left\|\pmb{\beta}\right\|^2, $$ którą można różniczkować względem parametrów $\pmb{\beta}$ i przyrównać do $0$. Dostajemy następujący estymator $$ \hat{\pmb{\beta}} = (\lambda I_N + X^\intercal X)^{-1}X^\intercal (Y - \hat{\beta}_0) $$
Iteracyjna - Gradient Descent
W tym podejściu do $S(\pmb{\beta})$ podchodzi się jak do funkcji kosztu nie wchodząc w jej interpretację i własności. Jest to problem z zakresu Convex Optimization.
Standaryzowania zmiennych dokonuje się poprzez odjęcie średniej i przeskalowanie do jednostkowej wariancji.
Dla próbki $x$ ustandaryzowana wartość $z$ obliczana jest w następujący sposób:
$$ z = \frac{x - u}{s},$$gdzie $u$ i $s$ to odpowiednio średnia i odchylenie standardowe populacji.
Dobrą praktyką jest standaryzować (ewentualnie skalować do przedziału) zmienne niezależne oraz zależne przed stosowaniem modeli. Dotychczas tego nie robiliśmy, ponieważ
W poprzednim zeszycie wspomniane zostało, że przy regularyzowanej regresji skalowanie jest niezbędne - skoro nakładamy taką samą karę na wszystkie parametry, to powinny one mieć wpływ w tej same jednostce na wszystkie zmienne niezależne.
Ponadto dla ustandaryzowanych zmiennych łatwiej jest porównywać modele między sobą, a także dobierać rozkłady prior dla parametrów, a tym będziemy się zajmować w tym zeszycie.
Należy jednak pamiętać, że współczynniki dla ustandaryzowanych zmiennych są cięższe do interpretacji. Dlatego wiele pakietów statystycznych domyślnie dokonuje standaryzacji i raportuje parametry dla ustandaryzowanych zmiennych, a także przeskalowane do oryginalnych rozkładów cech.
train_data.columns
W tym zeszycie intercept nie będzie traktowany, jak kolejna z cech. Nie chcemy go także standaryzować - otrzymalibyśmy wektor zer.
Pewien niepokój powinno budzi standaryzowanie zmiennych powstałych w ramach one-hotowania. Pomimo, że mają wartości liczbowe, to ciężko o nich myśleć jako o zmiennych ciągłych. Z drugiej strony chcemy zachować tą samą skalę dla wszystkich zmiennych w celu regularyzacji. Nie istnieje jednoznacza odpowiedź na to pytanie. Większość debaty w Internecie sprowadza się do
sklearn zwróci błąd jeżeli wywołamy
StandardScalerprzedOneHotEncoder
Należy pamiętać, że jesteśmy teraz poza światem założeń, więc najlepiej zewaluować oba podejścia w różnych warunkach i na tej podstawie wyciągnąć wnioski.
features = [
'Length3',
'Height',
'Width',
'species_0',
'species_1',
'species_2',
'species_3',
'species_4',
'species_5',
'species_6',
]
target = 'Weight'
scaler = StandardScaler()
train_data[features] = scaler.fit_transform(train_data[features])
test_data[features] = scaler.transform(test_data[features])
Zastosujemy implementację Ridge Regression z pakietu sklearn, która estymuje wyraz wolny w pkt 1. powyżej, czyli nie regularyzuje go. Posiada ona następujące solvery (dokumentacja):
‘auto’ chooses the solver automatically based on the type of data. ‘svd’ uses a Singular Value Decomposition of X to compute the Ridge coefficients. More stable for singular matrices than ‘cholesky’.
‘cholesky’ uses the standard scipy.linalg.solve function to obtain a closed-form solution.
‘sparse_cg’ uses the conjugate gradient solver as found in scipy.sparse.linalg.cg. As an iterative algorithm, this solver is more appropriate than ‘cholesky’ for large-scale data (possibility to set tol and max_iter).
‘lsqr’ uses the dedicated regularized least-squares routine scipy.sparse.linalg.lsqr. It is the fastest and uses an iterative procedure.
‘sag’ uses a Stochastic Average Gradient descent, and ‘saga’ uses its improved, unbiased version named SAGA. Both methods also use an iterative procedure, and are often faster than other solvers when both n_samples and n_features are large. Note that ‘sag’ and ‘saga’ fast convergence is only guaranteed on features with approximately the same scale. You can preprocess the data with a scaler from sklearn.preprocessing.
Wśród nich znajdziemy "analityczne" rozwiązanie w oparciu o operacje algebraiczne, a także interacyjne.
$\lambda$ nazywa się tutaj alpha :)
ridge_regression = linear_model.Ridge(
alpha=1.0,
fit_intercept=True,
normalize=False
)
ridge_regression.fit(train_data[features], train_data[target])
print(ridge_regression.intercept_) # nie bedzie sie zmieniac wraz ze zmiana alpha
pprint(dict(zip(features, ridge_regression.coef_))) # gdy damy nieracjonalnie duze alpha to beda same zera
ridge_regression.score(test_data[features], test_data[target])
(Zakładamy tutaj, że mamy już estymator wyrazu wolnego $\hat{\beta}_0$)
Przypomnijmy, że regresja liniowa to $$ (y_i - \hat{\beta}_0) = \pmb{x}_i \pmb{\beta} + \epsilon_i, $$ gdzie zakładamy $\epsilon_i \sim \mathcal{N}(0, \sigma^2)$. Czyli równoważnie możemy zapisać ten model jako: $$ (y_i - \hat{\beta}_0) | \pmb{x_i} \sim \mathcal{N}(\pmb{x}_i \pmb{\beta}, \sigma^2). $$
Log-likelihood takiego modelu to:
$$ \log p(D|\pmb{\beta}) = \frac{1}{N} \sum_{i=1}^N \log \left[ \frac{1}{\sqrt{2 \pi}\sigma}\exp \left( - \frac{((y_i - \hat{\beta}_0) - \pmb{x}_i \pmb{\beta})^2}{2 \sigma^2} \right)\right] = \text{const} - \frac{1}{2N\sigma^2}\sum_{i=1}^N \left| (y_i - \hat{\beta}_0) - \pmb{x}_i \pmb{\beta}) \right|^2 $$Gdy spojrzymy na Ridge Regression od strony MAP (nawiązanie do ostatnich zajęć) $$ \arg\max_{\pmb{\beta}} \log p(\pmb{\beta}| D) = \arg\max_{\pmb{\beta}} [ \log p(\pmb{\beta}) + \log p (D | \pmb{\beta})], $$ to okazuje się, że Ridge Regression to probabilistyczny model ze standardowym Gaussowskim priorem na wagach.
Niech $\pmb{\beta} \sim \mathcal{N}(\pmb{m}, V)$: $$ \log p(\pmb{\beta}) = - \frac{1}{2} (\pmb{\beta} - \pmb{m})^\intercal V^{-1} (\pmb{\beta} - \pmb{m}) + \text{constant},$$ gdy $\pmb{m}=0$ i $V = I_{m}$ to $$ \log p(\pmb{\beta}) = - \frac{1}{2} \|\pmb{\beta}\|_2^2 + \text{constant},$$ czyli otrzymaliśmy Ridge Regression!
Mamy zatem taki sam model, jak w przypadku Linear Regression, ale dla parametrów wybieramy rozkład prior. Na początek załóżmy, że wariancja błędu jest znana i skupmy się współczynnikach. Likelihood modelu wyrażony jest w następujący sposób $$ p(Y|X, \pmb{\beta}, \beta_0, \sigma^2) = \mathcal{N}(Y | \beta_0 + X \pmb{\beta}, \sigma^2 \mathrm{I}_N) \propto \exp \left(- \frac{1}{2\sigma^2} (Y - \beta_0 \pmb{1}_N - X \pmb{\beta})^\intercal (Y - \beta_0 \pmb{1}_N - X \pmb{\beta}) \right) ,\tag{1}$$ gdzie przez $\mathcal{N}(Y | \beta_0 + X \pmb{\beta}, \sigma^2 \mathrm{I}_N)$ oznaczamy prawdopobieństwem obserwacji $Y$ przy rozkładzie normalnym o parametrach $(\beta_0 + X \pmb{\beta}, \sigma^2 \mathrm{I}_N)$.
Na parametr $\beta_0$ (wyraz wolny) zakładamy nieinformatywny prior, czyli jego estymacja sprowadza się do MLE i tak jak wcześniej $$ \hat{\beta}_0 = \frac{1}{N} \sum_{i=1}^N y_i = \bar{Y}.$$ Dla ułatwienia notacji będziemy zakładać, że $Y$ zostało wycentrowane (wtedy możemy pominąć wyraz wolny), zatem $$ Y \leftarrow Y - \bar{Y}.$$
Wyśrodkujmy zmienną niezależną w naszym zbiorze danych.
y_mean = train_data[target].mean()
train_data[target] = train_data[target] - y_mean
test_data[target] = test_data[target] - y_mean
Zatem $(1)$ upraszacza się do $$ p(Y|X, \pmb{\beta}, \sigma^2) \propto \exp \left(- \frac{1}{2\sigma^2} (Y - X \pmb{\beta})^\intercal (Y - X \pmb{\beta}) \right)$$ które w dalszym ciągu jest Gaussianem, więc rozkład sprzężony to również Gaussian, oznaczmy go w następujący sposób $$ p(\pmb{\beta}) = \mathcal{N}(\pmb{\beta} | \pmb{\beta}_0, \pmb{V}_0)$$
Wtedy posterior obliczamy jako $$ p(\pmb{\beta}| X, Y, \sigma^2) = \mathcal{N}(\pmb{\beta} | \pmb{\beta}_0, \pmb{V}_0) \mathcal{N}(Y | X\pmb{\beta}_0, \sigma^2 \mathrm{I}_N) = \mathcal{N}(\pmb{\beta} | \pmb{\beta}_N, \pmb{V}_N),\tag{2}$$ gdzie $$ \pmb{\beta}_N = \pmb{V}_N \pmb{V}_0^{-1} \pmb{\beta}_0+\frac{1}{\sigma^2}\pmb{V}_N X^\intercal Y $$ $$ \pmb{V}_N = \sigma^2 (\sigma^2 \pmb{V}_0^{-1} + X^\intercal X)^{-1} $$ $$ \pmb{V}_N^{-1} = \pmb{V}_0^{-1} + \frac{1}{\sigma^2} X^\intercal X $$
Skąd to się wzięło? Murphy:eq4.125, Bishop:eq2.116
Ale my nie chcemy opuszczać "Bayesowskiej" drogi i chcemy doprowadzić do znalezienia rozkładu predykcyjnego (ang. predictive distribution). Możemy tak zrobić, ponieważ znamy pełen posterior, a nie tylko funkcję do niego wprost proporcjonalną.
Przez $\mathcal{D}=(X, Y)$ oznaczmy dotychczasowy zbiór treningowy, składający się z $N$ par. Chcemy znaleźć rozkład (predykcyjny) $y$ dla nowej próbki $\pmb{x}$. W ogólność chcemy wykonać następującą operację $$ p(y|\pmb{x}, \mathcal{D}) = \int p(y|\pmb{x}, \pmb{\beta}) p(\pmb{\beta}|\mathcal{D}) \mathrm{d} \pmb{\beta},$$ gdzie $p(\pmb{\beta}|\mathcal{D})$ to posterior (policzony na bazie priora i likelihoodu). W naszym wypadku jest dodatkowo sparametryzowany przez $\sigma^2$.
Zatem finalnie
$$ p(y|\pmb{x}, \mathcal{D}, \sigma^2) = \int \mathcal{N}(y | \pmb{x}\pmb{\beta}, \sigma^2) \mathcal{N}(\pmb{\beta} | \pmb{\beta}_N, \pmb{V}_N) \mathrm{d} \pmb{\beta} = \mathcal{N}\left(y | x\pmb{\beta}_N, \sigma^2_N(\pmb{x})\right),\tag{3}$$gdzie $$\sigma^2_N(\pmb{x}) = \sigma^2 + \pmb{x}\pmb{V}_N\pmb{x} $$
Warto zwrócić uwagę na $\sigma^2_N(\pmb{x})$, w której założona wariancja błędu jest powiększona o dodatkowy czynnik, który dodaje wariancji od posteriora wektora wag.
W przypadku probabilistycznego modelu nie mówimy już o przedziałach ufności i przedziałach predykcyjnych, tylko o credible intervals wynikających wprost z rozkładu predykcyjnego.
UWAGA 1: Taki model jest idealny dla uczenia przyrostowego. Dobry opis Bishop:p.154-156
UWAGA 2: Możemy znaleźć posterior dla modelu z priorem na $\sigma^2$ -> Murphy:p.236 rozdział 7.6.3 Bayesian inference when $\sigma^2$ is unknown
W powyższych rozważaniach mamy dwa hiper-parametry:
Oznaczmy $\eta= (\alpha, \lambda)$.
Popularnym rozwiązaniem tego problemu byłby Randomized Grid K-CV Search, gdzie dla losowego podzbioru konfiguracji z siatki przeprowadzona byłaby K-cross walidacja.
Alternatywą jest EB. Aby być w pełni Bayesowscy powinniśmy nałożyć hyperpriors na $\eta$. Wtedy możemy obliczyć posterior $\eta$ i rozkład predykcyjny ma następującą postać: $$ p(y|\pmb{x}, \mathcal{D}) = \int p(y|\pmb{x}, \pmb{\beta}) p(\pmb{\beta}|\mathcal{D}, \alpha, \lambda) p(\alpha, \lambda|\mathcal{D}) \mathrm{d} \pmb{\beta} \mathrm{d} \alpha \mathrm{d} \lambda$$
(Hiper-priory też mogą mieć swoje parametry i tak dalej...) Wycałkowanie $\alpha, \lambda$ jest analitycznie niemożliwe, obliczeniowo również. Pierwsze rozwiązanie to MAP - jeżeli posterior jest wyraźnie skupiony wokół wartości $\hat{\eta}$, to wybieramy ją i stosujemy jako plug-in estimate. By zaproponować alternatywę przypomnijmy jaką postać ma posterior $$ p(\alpha, \lambda | \mathcal{D}) \propto p(\mathcal{D}|\alpha, \lambda) p(\alpha, \lambda) $$ Jeżeli prior jest "płaski" i nie daje posteriora o skupionej masie, możemy zwrócić się ku MLE na likelihoodzie $p(\mathcal{D}|\alpha, \lambda)$. W obu przypadkach należy obliczyć posterior $$p(\mathcal{D}|\alpha, \lambda) = \int p(\mathcal{D}|\pmb{\beta}, \lambda) p(\pmb{\beta}|\alpha) \mathrm{d} \pmb{\beta},$$ a w przypadku MLE także przeprowadzić jego maksymalizację (opiera się na Expectation Maximization, które omówione zostanie w dalszej części semestru). Obie te rzeczy są dobrze opisane w Bishop:p.166-169 (UWAGA - inna notacja), rozdziały:
UWAGA: W tym przypadku online learning nie jest wprost aplikowalny - "empiryczna część" potrzebuje próbki by dokonać maksymalizacji (posteriora bądź evidence). Można spróbować Expectation Maximization w konfiguracji batchowej - ciekawy temat na projekt.
from matplotlib.animation import FuncAnimation
class UpdateBayesianLinearRegression:
'''Pomocnicza klasa do inkrementalnego uczenia'''
def __init__(
self,
ax,
model,
X_test_with_cat,
y_test,
x_idx,
x_label,
y_label,
category_idx=None,
interval=0.95
):
self.ax = ax
self.model = model
self.x_idx = x_idx
self.interval = interval
self.mean, = ax.plot([], [], color="red", label="Mean output")
self.credible_interval = ax.fill_between(
[],
[],
[],
color='cornflowerblue',
alpha=0.5,
label=f"Credible Interval"
)
self.num_samples = ax.text(0.5, 1, '', fontsize=15, horizontalalignment='right',
verticalalignment='bottom', transform=ax.transAxes)
self.ratio_in_interval = ax.text(1, 0, '', fontsize=15, horizontalalignment='right',
verticalalignment='bottom', transform=ax.transAxes)
if category_idx is not None:
for cat in np.unique(X_test_with_cat[:, category_idx]):
mask = X_test_with_cat[:, category_idx] == cat
self.ax.scatter(
X_test_with_cat[mask, x_idx], y_test[mask],
marker='+', s=100,
alpha=1,
label=cat
)
else:
self.ax.scatter(
X_test_with_cat[:, x_idx], y_test,
marker='+', s=100,
alpha=1,
)
# Set up plot parameters
self.ax.set_xlim(
np.min(X_test_with_cat[:, x_idx])-0.2,
np.max(X_test_with_cat[:, x_idx])+0.2
)
self.ax.set_xlabel(x_label)
self.ax.set_ylabel(y_label)
# Set up legend
handles, labels = self.ax.get_legend_handles_labels()
fig.legend(handles[:2], labels[:2], loc='upper right')
if category_idx is not None:
fig.legend(handles[2:], labels[2:], loc='lower right', title='Category')
def model_fit(self, i, X_train, y_train):
self.model.fit(X_train[i-1, :], y_train[i-1])
def __call__(self, i, X_train, y_train, X_test, y_test):
# This way the plot can continuously run and we just keep
# watching new realizations of the process
if i != 0:
self.model_fit(i, X_train, y_train)
predictive = self.model.predict(X_test)
predictive_mean = predictive.mean()
lower_interval, upper_interval = predictive.interval(self.interval)
ci_ratio = np.sum(
(lower_interval <= y_test) &
(y_test <= upper_interval)
) / len(y_test)
self.ax.set_ylim(
min(np.min(lower_interval), np.min(y_test))-0.5,
max(np.max(upper_interval), np.max(y_test))+0.5,
)
self.credible_interval.remove()
self.credible_interval = self.ax.fill_between(
X_test[:, self.x_idx],
lower_interval,
upper_interval,
color='cornflowerblue',
alpha=0.5,
label=f"Credible Interval"
)
self.mean.set_data(X_test[:, self.x_idx], predictive_mean)
self.num_samples.set_text(f'Number of train samples: {i}')
self.ratio_in_interval.set_text(f'Percent of observations within CI: {ci_ratio*100:.2f}%')
return self.mean, self.credible_interval, self.num_samples
W klasie poniżej zaimplementuj metodę predict tak by zwracała rozkład predykcyjny $Y$ dla zadanych danych wejściowych x. Wszystkie potrzebne wzory to $(2)$ i $(3)$, natomiast wszystkie potrzebne, nieznane parametry modelu zostały przypisane do zmiennych w konstruktorze oraz metodzie fit.
UWAGA: dodaj wyraz wolny do średniej
from scipy import stats
class ConjugateBayesLinReg:
def __init__(self, n_features, alpha, lmbda, intercept=0):
self.n_features = n_features
self.alpha = alpha
self.lmbda = lmbda
self.mean = np.zeros(n_features)
self.cov_inv = np.identity(n_features) * alpha
self.cov = np.linalg.inv(self.cov_inv)
self.intercept = intercept
def fit(self, x, y):
# If x and y are singletons, then we coerce them to a batch of length 1
x = np.atleast_2d(x)
y = np.atleast_1d(y)
# Update the inverse covariance matrix (Bishop eq. 3.51)
cov_inv = self.cov_inv + self.lmbda * x.T @ x # Vn-1
# Update the mean vector (Bishop eq. 3.50)
cov = np.linalg.inv(cov_inv) # Vn
mean = cov @ (self.cov_inv @ self.mean + self.lmbda * y @ x) #beta_n
self.cov_inv = cov_inv
self.cov = cov
self.mean = mean
def predict(self, x):
'''Metoda zwraca rozkład predykcyjny zmiennej niezależnej na bazie bieżących parametrów modelu.'''
# If x and y are singletons, then we coerce them to a batch of length 1
x = np.atleast_2d(x)
# TU WPISZ KOD !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
sigma = 1 / self.lmbda
y_pred_var = np.diagonal(sigma + x @ self.cov @ x.T)# sigma^2 + x Vn x
y_pred_mean = x @ self.mean + self.intercept
# Drop a dimension from the mean and variance in case x and y were singletons
y_pred_mean = np.squeeze(y_pred_mean)
y_pred_var = np.squeeze(y_pred_var)
return stats.norm(loc=y_pred_mean, scale=y_pred_var ** .5)
@property
def weights_dist(self):
return stats.multivariate_normal(mean=self.mean, cov=self.cov)
fig, ax = plt.subplots(figsize=(16, 8))
ud = UpdateBayesianLinearRegression(
ax,
model=ConjugateBayesLinReg(n_features=len(features), alpha=.3, lmbda=1),
X_test_with_cat=test_data[features+['Species']].to_numpy(),
x_idx=1,
x_label="Length3",
y_label="Weight",
category_idx=-1,
y_test=test_data[target],
)
anim = FuncAnimation(
fig,
ud,
fargs=(
train_data[features].to_numpy(),
train_data[target].to_numpy(),
test_data[features].sort_values(by=features[ud.x_idx], axis=0).to_numpy(),
test_data[[target, features[ud.x_idx]]].sort_values(by=features[ud.x_idx], axis=0)[target].to_numpy()
),
frames=train_data.shape[0],
interval=100,
blit=True,
)
plt.close()
anim